%% computeDataMoments.m
% This script computes the data counterparts of the moments appearing in
% Table 2. This script assumes the presence of a number of mat files with
% names "data_*" to be present in the working directory. These mat files
% have been added to the replication package and contain the data series as
% defined in Appendix C.1.1.

%% settings
startYear = 1970;
leverage = 1.5;

%% 1. capital stock data
capitalStockData = load('data_capital_stock');
startIdx = find(capitalStockData.year>=startYear,1,'first');
meanCapitalOutputRatio = mean(capitalStockData.capitalOutputRatio(startIdx:end));

%% 2. macro data
macroData = load('data_macro_aggregates');
startIdx = find(macroData.year>=startYear,1,'first');

[~, outputLogHP] = hpfilter(log(macroData.gdp(startIdx:end)),1600);
[~, consumptionLogHP] = hpfilter(log(macroData.consumption(startIdx:end)),1600);
[~, investmentLogHP] = hpfilter(log(macroData.investment(startIdx:end)),1600);

stdOutput = std(outputLogHP);
stdConsumption = std(consumptionLogHP);
stdInvestment = std(investmentLogHP);

meanConsumptionOutputRatio = mean(macroData.consumption(startIdx:end)./macroData.gdp(startIdx:end));
meanGovSpendingOutputRatio = mean(macroData.govSpending(startIdx:end)./macroData.gdp(startIdx:end));

corrConsumptionOutput = corr(consumptionLogHP,outputLogHP);
corrInvestmentOutput = corr(investmentLogHP,outputLogHP);

%% 3. surplus and debt data
surplusData = load('data_surpluses');
startIdx = find(surplusData.year>=startYear,1,'first');
startIdxLastDecade = find(surplusData.year>=2010,1,'first');
meanDebtOutputRatio = mean(surplusData.debtGdpRatio(startIdx:end));
meanDebtOutputRatioLastDecade = mean(surplusData.debtGdpRatio(startIdxLastDecade:end));
meanSurplusOutputRatio = mean(surplusData.surplusGdpRatio(startIdx:end));

[~, surplusOutputRatioHP] = hpfilter(surplusData.surplusGdpRatio(startIdx:end),1600);
[~, debtOutputRatioHP] = hpfilter(surplusData.debtGdpRatio(startIdx:end),1600);
stdSurplusOutputRatio = std(surplusOutputRatioHP);
stdDebtOutputRatio = std(debtOutputRatioHP);
corrSurplusOutputRatioOutput = corr(surplusOutputRatioHP,outputLogHP);

%% 4. asset price data
assetPriceData = load('data_asset_pricing');
startIdx = find(assetPriceData.year>=startYear,1,'first');
logYield5Year = log(1+assetPriceData.yield5Year(startIdx:end)/100); %transform percentage yield to log yield
bondReturnMonthly = 60/12*logYield5Year(1:end-1) - 59/12*logYield5Year(2:end); %approximate formula
equityReturnLeveragedMonthly = log(assetPriceData.equityIdx(startIdx+1:end))-log(assetPriceData.equityIdx(startIdx:end-1));

stdBondReturn = sqrt(12)*std(bondReturnMonthly);
stdEquityReturnLeveraged = sqrt(12)*std(equityReturnLeveragedMonthly);
meanBondReturn = 12*mean(bondReturnMonthly) + stdBondReturn^2/2;
meanEquityReturnLeveraged = 12*mean(equityReturnLeveragedMonthly) + stdEquityReturnLeveraged^2/2;
meanExcessEquityReturnLeveraged = meanEquityReturnLeveraged - meanBondReturn;
stdExcessEquityReturnLeveraged = sqrt(12)*std(equityReturnLeveragedMonthly-bondReturnMonthly);
sharpeRatio = meanExcessEquityReturnLeveraged/stdExcessEquityReturnLeveraged;

meanExcessEquityReturn = meanExcessEquityReturnLeveraged/leverage;
stdExcessEquityReturn = stdExcessEquityReturnLeveraged/leverage;

%% 5. risk-free rate data
riskFreeRateData = load('data_risk_free_rate');
startIdx = find(riskFreeRateData.year>=startYear,1,'first');
riskFreeRateMonthly = riskFreeRateData.riskFreeRate(startIdx:end);

meanRiskFreeRate = mean(riskFreeRateMonthly);
stdRiskFreeRate = std(riskFreeRateMonthly);

%% print results
fprintf('-------------------------------------------------\n')
fprintf('some moments in data:\n')
fprintf('-------------------------------------------------\n')
fprintf('                                                 \n')
fprintf('-------------------------------------------------\n')
fprintf('  std(Y):                     %.4f\n',stdOutput)
fprintf('  std(C)/std(Y):              %.4f\n',stdConsumption/stdOutput)
fprintf('  std(I)/std(Y):              %.4f\n',stdInvestment/stdOutput)
fprintf('  std(S/Y):                   %.4f\n',stdSurplusOutputRatio)
fprintf('  E[C/Y]:                     %.4f\n',meanConsumptionOutputRatio)
fprintf('  E[G/Y]:                     %.4f\n',meanGovSpendingOutputRatio)
fprintf('  E[S/Y]:                     %.4f\n',meanSurplusOutputRatio)
fprintf('  E[q^K K/Y]:                 %.4f\n',meanCapitalOutputRatio)
fprintf('  E[q^B K/Y]: full sample: %.4f, last decade: %.4f\n',meanDebtOutputRatio,meanDebtOutputRatioLastDecade)
fprintf('  E[r^E-r^B]:                 %.4f\n',meanExcessEquityReturn)
fprintf('  sharpe ratio:               %.4f\n',sharpeRatio)
fprintf('  corr(C,Y):                  %.4f\n',corrConsumptionOutput)
fprintf('  corr(I,Y):                  %.4f\n',corrInvestmentOutput)
fprintf('  corr(S/Y,Y):                %.4f\n',corrSurplusOutputRatioOutput)
fprintf('  std(q^B K/Y):               %.4f\n',stdDebtOutputRatio)
fprintf('  E[r^f]:                     %.4f\n',meanRiskFreeRate)
fprintf('  std(r^f):                   %.4f\n',stdRiskFreeRate)
fprintf('-------------------------------------------------\n')
